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ABSTRACT 

We analyze the link between active galactic nuclei (AGN) and mid-infrared flux using dust radiative 
transfer calculations of starbursts realized in hydrodynamical simulations. Focusing on the effect 
of galaxy dust, we evaluate diagnostics commonly used to disentangle AGN and star formation in 
ultraluminous infrared galaxies (ULIRGs). We examine these quantities as a function of time, viewing 
angle, dust model, AGN spectrum, and AGN strength in merger simulations meant to bracket the 
properties of ULIRGs. The more obscured burst begins SF-dominated with significant PAH emission, 
and ends with a ~ I0 9 yr period of red near-IR colors. At coalescence, when the AGN is most luminous, 
dust obscures the near-infrared AGN signature, reduces the relative emission from polycyclic aromatic 
hydrocarbons (PAHs), and enhances the 9.7/im absorption by silicate grains. Although generally 
consistent with previous interpretations, our results imply none of these indicators can unambiguously 
estimate the AGN luminosity fraction in all cases. Some identify relatively unobscured AGN where 
the direct torus emission is observed, while others indicate more highly obscured AGN. We show that 
a combination of the extinction feature at 9.7/um, the PAH strength, and a near-infrared slope can 
simultaneously constrain the AGN fraction and dust grain distribution for a wide range of obscuration. 
We find that this procedure, accessible to the James Webb Space Telescope, may estimate the AGN 
power as tightly as the hard X-ray flux alone, thereby providing a valuable future cross-check and 
constraint for large samples of distant ULIRGs. 

Subject headings: dust: extinction — Infrared: galaxies — Galaxies: interactions — Galaxies: star- 
burst — Radiative transfer 



1. INTRODUCTION 

Understanding the link between supermassive black 
holes (SMBHs) and their host galaxies is essen- 
tial for deciphering the formation and evolution of 
galaxies. Galaxy/SMBH co-evolution is expected 
given the observed correlation s between SMBH mass 
and galaxy properties (e-g-- IKormendv fc Richstonel 



1991 IMaeorrian et al l [19981 iFerrarese fe Merrittl I200C. 
Tremaine et al.H2002t : IGiiltekin et al.l 120091) and theoret- 



ical arguments that SMBH growth is self-regulated by 
feedback (e.g.. ISilk fc Reesl fl998t ISpringel et all [20051: 
IHopkins et alJl2007allblL " 

A key prediction from this framework is that a 
galaxy experiencing rapid inflow of cold gas can evolve 
through various classes of starbursts and active galac- 
tic nuclei (AGN), such as ultra- luminous infrared galax- 
ies (ULIRGs) and quasars (QSOs), and that these 
phases are connecte d in an evolutionary sequence (e.g., 
I Sanders et al.l I1988D . With such a model, both a 
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starburst and a he avily obscured AGN can co-exist 
([Hopkins et al.ll2006l ). 

Testing this picture requires not only finding signatures 
of obscured AGN activity within starburst galaxies, but 
also interpreting them in the context of galaxy/ SMBH 
co-evolution. Locally, the most extreme starbu rsts are 
the ULIRGs (L m > 10 12 £^ : [Sanders et al.lfl988h . which 
are almost e xclusively the result of a rec ent / ongoing ma- 
jor merger ([Sanders fc Mirabel Il996af ). which triggers 
both starburst and obscured AGN activity. To evaluate 
SMBH growth and the role of feedback during this crit- 
ical phase, it is necessary to robustly estimate the AGN 
power during all observed phases. A primary challenge to 
determine the fraction of flux attributable to the AGN 
is the reprocessing of its photons by the host galaxy's 
interstellar medium (ISM). In essence, most diagnostics 
focus on finding signatures that are directly associated 
with the AGN output: narrow- line region (NLR) emis- 
sion, X-rays, torus emission, or radio emission. 

A popular approach utilizes the mid-infrared (mid- 
IR) regime, loosely defined as 3-30/xm, where typical 
starbursts have a spectral energy distribution (SED) 
with a minimum in continuum emission, while contin- 
uum emission from dust surrounding AGN rises toward 
longer wavelengths owi ng to hot dust emission, poten- 
tially from a torus (e.g. , 
[2^[H5niTet^[2006l: 

This idea led to the construction of diagnos- 
tic diagrams separating starburst from AGN domi- 
nated sources, based on Infrared Space Observatory 
(ISO ) data of local IR-luminqus sources dLutz et all 
[19961 [19981 IGenzel et all H998I: IRigopoulou et all 119991: 
lLaurent et al.l 120001: iTran et al.l 120011) . The sensitive 
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Infrared Spectrograph (IRS; IHouck et al . 2004) aboard 
the Spitzer Space Telescope ( Werner et al.l 1 2 0041) en- 
abled the expansion of this approach to larger sam- 
ples covering a wider range of physica l properties (e.g.. 
Smith et al.l 120071: IBrandl et al l 120061: ISchweitzer et alj 
2006jlWu et al.ll2006l:IDale et aU2006|:IArmus et al.ll2006l 
20071: iSturm et alj|2006t iSpoon et aTll2007f ) andpermit- 
ted mid-IR diagnostic work on objects at higher red- 
shifts. In particular, mid-IR diagno stics have been ap- 
plied to g alaxies at z ~ 1 -3 (e. g. IHouck et al.l 120051 : 
lYan et all 120071 : iPope et al.1 12008T) the epoch of peak 
star- formation rate density (e.g., iBouwens et al.l 120071) 
and p eak quasar number density (e.g.. iRichards et al.l 
l2006bf) . which makes it a period crucial for evaluat- 
ing the co-evolution of galaxi es and SMBHs. It is also 
when LIRGs and ULIRGs (see lSanders fc Mirabellll996bL 
for a review) make a significant contribution to the 
global averaged lumi nosity density and to the cosmic 
infrared background (iLe Floc'h et al.l 120051: iDole et al.l 
l200llCaputi et al.ll2007t IHopkins et al.ll2010D . 

Other findings on high-z IR-luminous galaxies suggest 
that high-redshift ULIRGs may be unlike local ULIRGs 
in that they may not all be late-stage mergers. For exam- 
ple, z ~ 2 sources tend to be more starburst-dominated 
than z ~ sourc es of comparable mid-IR luminosity 
(|Fadda et al.ll2010T ). Differences in both the spectral en- 
ergy distributions (SEDs) and morphologies of high red- 
shift ULIRGs, as compared with local (z < 0.1) sources, 
do su ggest that they are not analogous fe.g.. lSaiina et al.1 
12012L hereafter S12). Furthermore, submillimetre galax- 
ies, a subset of high-redshift ULIRGs, may be a mix of 
early-stage quiescently star-forming mergers, late-stage 
merger-induced starbu rsts, and isolated disk galaxies 
(lHavward et alJl2012alfbh. The typically greater gas frac- 
tions (|Tacconi et al.ll2010| ) imply that fueling both star- 
formation and black hole accretion is relatively easier, 
and therefore not necessarily analogous to what we see 
in local ULIRGs. 

In both the local and high-redshift work, it is not clear 
to what extent common diagnostics miss the crucial evo- 
lutionary stage when the AGN is most deeply obscured. 
It is possible that we can only detect the presence of 
an AGN once sufficiently many unobscured lines of sight 
to the AGN exist. Moreover, while IR diagnostics gen- 
erally follow the expected AGN/starburst fractions, the 
evolution of IR SED properties reflects a c omplex mix 
of these components (jYounger et al.l l2009bf l which may 
be di fficult to interpret f or individual sources or sam- 
ples (jVeilleux et al.l 120091 hereafter V09). Signatures 
thought to be evidence for AGN, such as a power-law 
mid-IR SED, may also b e caused by intense starbursts 
([Narayanan et all l2010bl ) . Therefore it is desirable to 
search for signposts of AGN powering that may suffer 
minimally from these complexities. 

In this paper we constrain the applicability of mid-IR 
diagnostics and test their common interpretations, in- 
tending to gain insight regarding AGN activity inferred 
by conventional methods. In Section [2] we describe our 
3D dust radiative transfer calculations of two representa- 
tive starbursts meant to bracket the plausible conditions 
in ULIRGs, from which we compute and analyze fea- 
tures of the mid-IR SEDs as described in Section [3J In 
Section [4[ we compare these diagnostics to observations 
and evaluate how they depend on AGN power, evolu- 



tionary stage, viewing direction, intrinsic AGN SED, and 
assumptions about the dust and ISM. We consider im- 
plications of this work for future studies of the evolution 
of IR-luminous galaxies in Section |6l and we conclude in 
Section [7] 

2. SIMULATIONS 



We combine high-resolution Gadget-2 (|Springell 
2005) 3-D N-body/smoothed-particle hydrodynamics 
(SPH) si mulations of equal-mass galaxy merge rs with the 
Sunrise (| Jonssonl [20061 : Uonsson et al.1 l2010f ) polychro- 
matic Monte Carlo dust radiative transfer (RT) code to 
examine commonly-applied mid-infrared (mid-IR) AGN 
signatures. We focus on two representative mergers, a 
"highly obscured" hyper-LIRG analogue, and a "less ob- 
scured" marginal ULIRG example. The highly obscured 
simulation presented here is the same one focused on in 
lHavward et al.l ([201 ll ). so here we summarize the essen- 
tial assumptions and differences from those works. 

2.1. Hydrodynamical simulations 
Gadget-2 is a TreeSPH ([Hernquist fc Katzl 



1989) code that conserves both energy and entropy 
( Springel fc Hernquistl l2002f) . The simulations i nclud e 
radiative heating and cooling as in IKatz et all (fl996h . 
Star formation (SF) is m odeled to reprod uce the global 
Kennicutt-Schmidt law ([Kennicuttl 119981 ) via the vol- 
umetric relation psfr oc Pgks with a minimum density 
threshold n ~ 0.1 cm~ 3 . We do not track the formation 
of molecular gas or resolve individual molecular clouds. 
Thus the KS law employed should be considered simply 
an empirically and physically motivated prescription to 
summarize physics we do not resolve. However, this 
and other details of the hydrodynamical simulations 
are relatively unimportant for our analysis; all that we 
require are galaxy models that approximately resemble 
real galaxies on which we can perform radiative transfer. 

The structure of the ISM is modeled via a two- 
phase sub-resolution model in which cold, star-forming 
clouds are embedded in a diffuse, hot medium 
(|Springel fc Hernquistl |2003|) pressurized by supernova 
feedback th at heats the diffus e ISM and evaporates the 
cold clouds ([Cox et al.ll2006bl ). Metal enrichment is cal- 
culated by assuming each gas particle behaves as a closed 
box with a yield y = 0.02. Supermassive black hole 
(SMBH) particles accrete via Eddington-limited Bondi- 
Hoyle accretion and deposit 5% o f their luminosity to 
the nearby ISM as the rmal energy ([Snringel et al.l [20051 : 
iDi Matteo et al.l 120051 ). The luminosity is computed 
from the accretion rate assuming 10% radiative efficiency, 

Lbol = O.lAfc 2 . 
Recent work by iVogelsberger et all ([2012T ). 



Siiacki et al 



Torrev et al 



iKeres et al. 
has shown that 



Gadget 



and 
can, 



in some circumstances, yield results that disagree with 
simulations performed with the moving mesh code 
Arepo (Springel 2010). These comparisons involved 
cosmological simulations of galaxy formation, where 
the phase structure of the gas can be complex on small 
scales. The isolated merger simulations described here 
involve gas that has simpler physical properties because 
the sub-resolution model we employ prevents it from 
developing multiple phases in the ISM locally. For 
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Figure 1. An overview of the SEDs of our merger simulations. On the left we show false-color rest-frame U-V-J and IRAC composite 
images of the highly obscured merger at each of four example times, as seen from the same direction as the SEDs to the right are 
measured, and utilizing our default RT model and AGN strength. In the middle two columns we show total rest-frame SEDs of the highly 
obscured merger from an example viewing direction at each of four times labeled in the left-most panel. SEDs have the same arbitrary 
flux normalization. Different curves in the left column represent our fiducial dust model, but the input AGN spectrum has been multiplied 
by 0, 1, and 10. Curves in the second column show our fiducial AGN level, but different ISM assumptions. In the right two columns we 
similarly demonstrate the SEDs of the less obscured merger. 

this reason, explicit comparisons between Gadget 
and Arepo in the context of merger simulations like 
those we study here (C. Hayward et al., in preparation) 
demonstrate that the conclusions discussed in this paper 
are not affected by the choice of simulation method. 



2.2. Radiative transfer 

We use Sunrise in post-processing to calculate the 
SED observed from seven cameras distributed isotrop- 
ically in solid angle every 10 Myr. For a full descrip- 
tion o f Sunrise see Uonssonl (|2006l ) and Uonsson et al.l 
d2OT0h . and for a summary of components essential for 
this work and a descr iption of other specific choices, see 
lHavward et "all ff20TTh . 

Sunrise calculates the emission from the stars and 
AGN in the Gadget- 2 simulations and the atten- 
uation and re-emiss ion from dust. Starburst99 
(jLeitherer et al.|[l999h SEDs are assigned to star parti- 
cles according to their ages and metallicities, and SMBH 
particles emit the lu minosity-dependent templates of 
Hopk ms et~aLl ([20023) by assuming the formula above 
for I/boi- At mid-IR wavelengths this template emits the 
mean SED of lRichards et alj (|2006aD . In Section SH we 
vary the AGN source SEP, ap plying two clumpy torus 
models bv lNenkova et all i|2l)l)Si that span their modeled 
properties. 

To calculate the dust density, Sunrise projects the 
Gadget-2 gas-phase metal density onto a 3-D adaptive 
mesh refinement gr id. Here we have assumed 40% of the 
metals are in dust (Dwck 1998). We consider the Milky 
Way R v = 3.1 (MW) and SMC bar (SMC) dust models 
of lWeingartner fc Draind (|2001l) updated bv lDraine fc Lil 



((20071) . 

Sunrise then performs Monte Carlo radiative transfer 
by emitting photon packets from the sources and drawing 
interaction optical depths from the appropriate probabil- 
ity distribution as the packets traverse the ISM. For each 
grid cell, the temperature of each dust species (with the 
exception of polycyclic aromatic hydrocarbons; PAHs) is 
calculated assuming the dust is in thermal equilibrium, 
and the dust re-emits the absorbed energy as a modified 
blackbody. Half the luminosity absorbed by PAHs with 
size a < 100 A is reemitted assuming the grains are in 
thermal equilibrium. The other half of the luminosity 
is emitted as a fixed observati onally-derived temp late (a 
sum of Lorentzian profiles; see lJonsson et al.1l2010l for de- 
tails). In high-density regions, the dust can be opaque to 
its own emission, so the contribution of the dust emission 
to dust heating must be considered. Sunrise computes 
this self-consistently by iteratively performing the trans- 
fer of the dust emission and the temperature calculation 
using a reference field technique. 

Our focus herein is the mid-IR portion of the rest-frame 
SEDs, 1 fun < A < 30/im, a regime in which dust grains 
heated to a wide range of temperatures (T ~ 100-1500 
K) emit thermal radiation. Arbitrary dust tempera- 
ture distributions are possible during a gas-rich galaxy 
merger, owing to the complex and chaotically-evolving 
multi-phase structure of the ISM induced during a global 
starburst and feedback-regulated SF and SMBH growth. 
Thus the accurate dust heating calculations employed by 
Sunrise, including treatments of multiple dust species, 
are essential for understanding this key phase of galaxy 
evolution. Moreover, studies such as this one are crucial 
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for verifying these methods. 

Since we will make extensive use of the 9.7 /jm sili- 
cate absorption feature, we stress that the depth of the 
feature is determined self-consistently via the radiative 
transfer. The strength of the feature depends on the 
amount of dust and the geometry of sources (stars and 
AGN) and dust. For reasonable galaxy geometries we 
expect the AGN to be more obscured than the bulk of 
the stars; this differential extinction is inherently ignored 
when one assumes a foreground screen dust geometry, so 
our simulations provide an excellent way to test the un- 
certainties caused by this assumption. 

The results of the Sunrise calculation are spa- 
tially resolved, multi-wavelength SEDs observed from 
seven directions. The success of this approach at 
modeling dive r se gala x y populations — bo t h local (e.g. . 
Younger et all I2009bt [Bush et all l2010t ISnvder et al 



201 II) and hig h -redshi ft (e.g.. IWuvts et all 1201 



Narayanan et 



high -redshi rt (e.g.. IWuvts et al.l IzUJ 
al] 12010113 IHavward et al.l !2012alM 



Si 



lends credibility to its application here. 

2.3. Galaxy models 

We focus here on two mergers meant to bracket the 
range of plausible conditions in ULIRGs. 

In each case, we use two identical progenitor galax- 
ies. The progenitors of the highly obscured merger are 
composed of an exponential disk with baryonic mass 
4 x lO n M (60% of which is gas) and dark matter halo 
of mass 9 x 1O 12 M described by a iHernquistl (|1990f ) 
pro file. The galaxy propert ies are scaled to z ~ 3 follow- 
ing [Robertion|eF^y ()2006D . The or bital param e ters ar e 
identical to those of the 'e' orbit of ICox et al.l ([2006a). 
We refer the reader to fu rther details of the merger setup 
in IHavward et aLl (|2011fl . This merger scenario is unique 
only insofar as it induces a massive, gas-rich starburst 
with elevated SMBH accretion. 

This highly obscured merger reaches a peak log Li r ~ 
13 briefly, having \ogLm > 12 for ~ 1 Gyr. In con- 
trast, the progenitors of the weakly obscured example 
have lower mass (M» ~ 1O 1O M ) and are less gas-rich 
(40%) than the highly obscured example. This example 
briefly reaches Lm ~ 1O 12 L under our default assump- 
tions. Our discussion will focus mainly on signatures 
of AGN activity that normalize out the galaxy's mass 
or total luminosity, and therefore the salient difference 
between these two cases are the gas-richness and (conse- 
quently) the relative extents to which the merger triggers 
star formation, AGN activity, and obscuration. 

These "high obscuration" and "low obscuration" cases 
do not necessarily accurately reflect the real ULIRG pop- 
ulation at a given redshift. The two experience very dif- 
ferent physical conditions during their starbursts, and 
for any given simulation of this type, the physical state, 
in particular the ISM structure, is rather uncertain. In 
some ways, this is realistic, as the population of ULIRGs 
may be a heterogeneous one, with contributions from 
sources of different types. The "high obscuration" ex- 
ample is, by choice, a particularly massive and gas-rich 
merger, and thus becomes a particularly luminous hyper- 
LIRG. While this kind of source occurs at higher redshift, 
it is uncommon in the local universe. In contrast, the 
"weakly obscured" example is marginally a ULIRG, and 
may be very common at all redshifts. 



In Figure Q] we summarize the observed SEDs of our 
mergers. We focus on four times during each merger 
spanning ~ 500 Myr, and present the total mid-IR SEDs 
emerging in a particular direction. In the left SED 
columns we show the results of our fiducial ISM assump- 
tion, but artificially vary the total luminosity of the two 
(and eventually one) SMBH particles. In the right SED 
columns we explore our assumptions about the ISM and 
dust, a study that we analyze in more detail below in 
Section [O] (see also Section 12.41) . 

Alongside the SEDs we show high-resolution false- 
color composite images corresponding to the same view- 
ing direction as the highly obscured SEDs. The left 
images present the system in the U, V, and J bands 
(| Johnson fe Morganlll953T) . and the second in channels 
one, tw o, and four of the Spitzer Infrared Array Camera 
(IRAC; iFazio et al.l 12004) . These images provide some 
insight into how these systems might be classified in 
terms of a merger-stage diagnostic - if this system is ob- 
served at high redshift, then the only time it is obviously 
a merger is at t < t\. This precedes the SMBH's peak 
luminosity and therefore the times at which it makes an 
obvious contribution to the mid-IR SED, a period span- 
ning ~ 1 Gyr. 

2.4. Alternate Radiative Transfer Models 

In order to gain additional physical insight, we con- 
sider different sets of assumptions for the radiative trans- 
fer stage, including two sub-resolution models for the 
dust distribution on scales below that resolved by the 
Gadget- 2 simulation. For our highly obscured simula- 
tion, we will focus on the ISM treatment referred to in 
IHavward et al.l (|201 If ) as "multi-phase off" . This choice 
assum es that the dust mass conta ined in both phases 
of the ISpringel fc Hernquistl (|2003| ) ISM model is dis- 
tributed uniformly across each resolution element, which 
is likely appropriate for gas-rich mergers in which the 
central regions are are composed of dense, almost exclu- 
sively molecular gas. Throughout this work we refer to 
this model as our "default ISM" treatment. 

Alternatively, we can assume that gas in th e cold, 
dense clouds of the ISpringel fc Hernquist] (|2003l ) model 
has negligible volume filling factor. This alternate choice 
distributes dust mass that corresponds to the diffuse 
gas phase of the ISpringel fc Hernquistl (|2003t ) model uni- 
formly across each resolution element, but assumes no at- 
tenuation from the dust occupying the dense, cold clouds. 
In this case we neglect attenuation and emission from 
both the cold clouds in which stars are formed and from 
other clouds photons encounter along the line-of-sight. 
This assumption thus gives a lower limit on the amount 
of attenuation. We refer to this assumption as our "al- 
ternate ISM" treatment, and apply it to both our highly 
and weakly obscured simulations. 

The two treatments of subresolution dust structure 
we employ should be considered two plausible and 
physically-motivated, yet uncertain, ISM models that 
encapsulate unresolved processes. With current simula- 
tions it is only possible to parameterize our ignorance in 
this way. However, any model used to interpret ULIRG 
SEDs is limited by this uncertainty. These simulations 
are thus particularly useful for their ability to quantify 
the uncertainty caused by dust dumpiness. 

For both mergers, we perform the RT calculations af- 
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ter multiplying the SMBH luminosity at all times by zero 
(AGNxO), one (AGNxl), and ten (AGNxlO), allowing us 
to manually adjust the AGN contribution to the SEDs 
(however, the effects of AGN feedback are kept fixed). 
This enables us to cleanly test how a given indicator 
depends on AGN luminosity for fixed geometry. The 
AGN contribution to the mid-IR SED arises from sepa- 
rate self-consistent RT calculations with each of our as- 
signed AGN luminosities, dust grain models, and ISM 
assumptions. At each of these AGN strengths, we apply 
both ISM assumptions to the highly obscured merger and 
the "alternate ISM" treatment for the less obscured case. 
For these three sets of three simulations, we perform the 
RT using both the MW and SMC bar dust models. 

2.5. X-ray Calculations 

We compute X-ray fluxes from SMBH particles 
and subsequent attenuatio n by the ISM foll owing 
the approach describ ed by iHopkins et al.l (|2005j ) and 
Hopkins et al.l (|2006f ) . This m ethod uses the intrin - 
sic quasar SED ( Section 12.21 and [Hopkins et al.l 12007a ). 
and calculates its extinction at X-ray wavelengths by 
applying the photoelectr ic abso rption cross sections of 
iMorrison fc McCammonl (|1983f ). as well as Compton 
scattering cross sections, eac h scaled by met a llicity . 

A key difference between Hopkin s" et al.l (|2005) and 
the present work is their assumptions correspond to the 
"alternate ISM" treatment we described in Section 12 .41 
where the cold clumps in the ISM, and hence often a sig- 
nificant fraction of the dust mass, are assumed to have a 
small volume filling factor. This may be true under many 
conditions, but may not be applicable to extremely gas- 
rich ULIRGs with abundant supplies of molecular gas. 
Therefore for our highly obscured merger simulation we 
calculate the column densities that attenuate the X-ray 
flux both ways: by discarding the cold phase mass ("al- 
ternate ISM"), and by keeping the cold phase mass ("de- 
fault ISM" ) . For our less obscured merger simulation, we 
use only the "alternate ISM" model, discarding the cold 
phase mass. 

2.6. AGN Fraction 

For this paper, we focus on the AGN fraction, defined 
to be the total power emitted by the SMBH particles 
divided by the total power emitted by all sources in 
the simulation, over the wavelength range used by Sun- 
rise: 0.09-10Vm. We denote this by L agn /£boi- This 
quantity does not include power emitted at X-ray wave- 
lengths, but we have checked that all subsequent analysis 
does not depend on this restriction, owing to the tight 
correlation between X-ray and non- X-ray luminosity of 
the input AGN source. 

This quantity is not normally available for a given 
observed sample. Here, the SMBH accretion rate and 
-^agn/^boi are available directly from the Gadget and 
Sunrise calculations. Thus we can evaluate directly 
the effectiveness of observed indicators, defined in Sec- 
tion [3j at estimating the AGN contribution utilizing self- 
consistent calculations of the galaxy's stars, dust, and 
SMBH. 

3. MID-INFRARED SPECTRAL DIAGNOSTICS 
3.1. From Simulations 




X (|i m) 

Figure 2. In the solid black line we show an example SED out- 
put from Sunrise in the source's rest frame. For all subsequent 
analysis we use this low-resolution version, but here we plot in 
blue a higher-resolution calculation that demonstrates the validity 
of the continuum estimation that we show in gray (see Section (3{. 
In the dashed black curve we show the intrinsic AGN spectrum 
used by Sunrise as a light source from the SMBH particles. This 
SED occurs between the peak of the starburst and the peak of the 
AGN's bolomctric contribution to the SED (between times t2 and 
£3 shown in Figure [TJ . Note that during this time the intrinsic 
AGN emission (attributed to the torus) at 2fim < A < lOfim is 
completely absorbed and reprocessed into the far-infrared. 

We compute mid-infrared diagnostics based on the 
rest-frame spatially integrated Sunrise SEDs, denoted 
juW) or /^(A). For simplicity and to follow common 
practice, we denote /„(A) by fx/^m- Owing to compu- 
tational constraints, these SEDs have very low spectral 
resolution (A/AA ~ 10), so we limit consideration to 
approximations of several observationally-motivated es- 
timators. We focus on two spectral dust features, rg.7, 
ry.r, defined by 

/9.7 

rg.7 = 

eg. 7 



where c\/^ m is the continuum level estimated by inter- 
polating linearly between fejj and /15. The remaining 
indicators are ratios of f\ at 1.6, 6, 15, and 30 nm. We 
present all SEDs and SED quantities in the source's rest 
frame. 

In Figure [5] we present an example rest- frame SED in 
the mid-IR and highlight the spectral features that we 
will consider here. In addition to the final SED, we plot 
the intrinsic SMBH SED that we assume, which exhibits 
a generally flat shape with an "IR bump" . If this source 
were observed unobscured, or obscured only by an opti- 
cally thin (at several /im) dust column, then this would 
resemble a "power-law-like" AGN source with red near- 
and mid-IR colors. We see this situation in Figured) the 
less obscured merger has drastically different /4.5//1.6 
slopes depending on the AGN strength assumed, and a 
comparison of the "no dust" and "default ISM" cases 
indicates that dust attenuation is insignificant at wave- 
lengths longer than a few microns. Furthermore, this 
property of AGN SEDs has be en used to select AGN i n 
wide-field IR surveys (see, e.g., IStern et al.l [20051 12011) . 
Thus we will keep this feature in mind while analyzing 
the mid-IR properties of more highly obscured sources. 
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By contrast, the starburst SED falls as lambda in- 
creases at lfim < A < 10 fim, but also exhibits dis- 
tinctive features from dust grains identified as polycyclic 
aromatic hydrocarbons (PAHs) superimposed at 3.3, 6.2, 
7.7, 8.6, 11.3, 12.6, and 17/xm. 

3.2. From Observations 

Obscrvationally, weaker PAH and stronger continuum 
(i.e., lower PAH equivalent width (EW)) effectively in- 
dicate a higher relative contribution fro m the AGN in 
dusty galaxies, especia l ly ULIRGs (e.g. , ILaurent et al.l 
[^OOaiSturm et al.|[2lMlTran et al.ll200l . This diagnos- 
tic is supported by mid-IR fine structure line diagnostics 
implying a n AGN-like radiation field in sources o f lower 
PAH EW (IGenzel et al.iri998t lArmus et al.l[2007h . 

More recently, such mid-IR diagnostic studies dis- 
covered sources that have a small PAH EW but a 
deep 9.7/xm silicate absorption feature (e.g.. lHouck et al.1 
12005ft . which seems to req uire a deeply embedd ed, cen- 
trally concentrated source (|Levenson et al.ll2006f) . These 
sources are believe d to represent Compton - thick AGN 
(|Bauer et al.l 120101 : iGeorgantopoulos et al.l 1201 ID . al- 
though there is at least one known star forming dwarf 
without an AGN whose mid -IR spectrum fits the above 
criteria (jRoussel et al.|[200l . 

In Section [SJ we will compare our predictions with 
data for local and high-z starbursts with available low- 
resolution Spitzer IRS spectra. Specifically, we use the 
ULIRG sample of V09 and a sample of 192 24/im-selected 
z^0.3-3 starbursts and obscured quasars. The red- 
shifts and I RS spectra of the 24/ zw-bright sample ar e 
presented in lYan et al.l (|2007f ) and iDasvra et all (120091) . 
while the full IR SEDs are compiled in iSaiina et al.l 
(HoH) (S12). 

Using the IRS spectra for both ULIRG samples, we 
estimate rr.r, rg.7, fe, /15, and /30 in the same manner as 
for the simulated SEDs. However, the observed spectra 
have a more limited spectral coverage. The observed 
coverage of the IRS spectra varies between 5.2 and 38/irn, 
or 14-38yLtm if only the LL modules are used, as is the 
case for most of the higher- z sample. This means that 
for the local ULIRG sample, the required rest-frame 6- 
30/zw is covered by the available IRS spectra, while for 
the 24/im-bright sample, the 15/um and 30/jm continuum 
points are outside the IRS coverage above z ~ 1.5. Below 
z ~ 1.5, the 6/im continuum point is often outside the IRS 
coverage (depending on whether or not the SL module is 
available), and even the 7.7/jm point can be outside the 
IRS coverage below z~0.8. We compensate us ing the 
empirical SED model fits from lSaiina et all (|2012f) , which 
necessarily introduces additional systematic uncertainty. 
We constrain these extrapolations below or above the 
IRS spectral coverage with IRAC 3.6, 4.5, 5.8, and 8/j,m, 
and MIPS 70/im photometric points, ensuring that the 
mid-IR colors obtained are reasonably accurate. 

We choose to exclude the z < 0.9 sources from this sam- 
ple in order to avoid all cases where the 7.7/jm PAH fea- 
ture is not covered by the IRS spectra. This step also re- 
moves all sources whose overall IR luminosities (£3-1000) 
place them in the LIRG rather than ULIRG category. We 
also exclude the z > 2.5 sources, for which the 9.7 /jm sil- 
icate absorption feature is too poorly defined, being only 
partially covered by the IRS spectra. This defines a total 
of 118 sources. We explore determining the above quanti- 
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Figure 3. We compare various rest-frame mid-IR quantities as a 
function of time for our highly obscured (left) and less obscured 
(right) starburst simulations. At each time, we plot the median 
over viewing angle of each quantity. Vertical lines highlight the 
four characteristic times t\-t^ from Figure[TJ corresponding to pre- 
burst, peak starburst, peak AGN, and post-burst. The three curves 
in each panel represent the same models shown in the left column 
of Figure [TJ we multiply the AGN luminosity deduced from the 
simulation by 0, 1, or 10. In several rows, the solid gray (patterned 
blue) boxes subtend the 25-75 % range in the distri bution of local 
ULIRGs (quasars) compiled by Vcilleux et al. (2009). In the highly 
obscured case, the AGN is most dominant and most obscured at 
the time £3, when the PAH emission (^7.7) and Si absorption (rg.7) 
ratios are minimized, in agreement with previous reasoning. The 
mid-infrared colors /30//15 and fis/fe also exhibit a weaker sig- 
nal around this time, but this signal can be confused with the 
starburst phase (£2) for fis/fe, and depends non-linearly or not at 
all on the AGN strength with /3o//l5- The more typical starburst 
represented in the right column is a marginal ULIRG, and its mid- 
infrared SED is relatively insensitive to the presence of the AGN 
which dominates at 1.4 < t/Gyr < 1.8. 



ties with and without interpolating the IRS spectra onto 
the much lower resolution simulated SED wavelength ar- 
ray, and we find no significant difference in the overall 
results. 

Lastly, in 4 cases, the silicate feature is saturated mak- 
ing the measured 7*9.7 value a lower limit. The saturation 
flag is triggered when the mean of the IRS spectra be- 
tween 9.0 and 10.4/im is less than or equal to the stan- 
dard deviation in the same region. 

4. SIMULATION RESULTS 
4.1. Time Evolution 
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In Figure |31 we explore how the mid-IR diagnostics 
vary with time for our two fiducial simulations: "highly 
obscured" and "less obscured". We show spectral diag- 
nostics along with quantities of interest such as the over- 
all IR luminosity, the star formation rate (SFR), and the 
AGN fraction, L agn /Lho\- We assume our default ISM 
model and default AGN torus model, but we show the 
effect of varying these choices in Section 14.21 and Sec- 
tion !4.4[ respectively. All curves shown are the median 
over the viewing angle. The scatter introduced by the un- 
certainty in the viewing angle is addressed in Section [4.3l 

To guide the discussion, we indicate the four times from 
the panels in Figure [T] with vertical lines for the highly 
obscured (left column) and less obscured (right column) 
mergers. The progenitors are still separated at time t\, 
which precedes the period of final coalescence and peak 
Lir by ~ 200 Myr. The maximum SFR and AGN power 
occur at t ~ t^-t-i, after the galaxies have merged, but 
the central sources are still largely obscured by dust. At 
i4, up to ~ 300 Myr later, the SFR has plummeted but 
the AGN fraction is still elevated, and the remnants are 
on their way to being ellipticals. 

Qualitatively, we find that our simulated mid-IR diag- 
nostics behave largely as expected, but not always. Be- 
low we address each of the diagnostics in turn. 

The PAH strength, r?.7, drops at the time of coa- 
lescence fa—tz), reaching its lowest values at £3 when 
the AGN fraction is highest. Pre-coalescence, the PAH 
strength is weakest for the model with strongest AGN 
(AGNxlO). Post-coalescence, the PAH strength again 
increases but does not return to its pre-coalescence lev- 
els. No significant "coalescence" dip in the PAH strength 
is seen in the less obscured case. The level of absorption 
by silicate grains, rg.7, is greatest when the AGN is most 
luminous, which is also when the relative strength of the 
PAH emission is weakest. These extrema are relatively 
sharp and confined to the period of final coalescence at 
~ ti-tz, although the AGN luminosity fraction remains 
elevated for w 1 Gyr. Overall, the shape of the PAH 
strength curve for either the highly obscured or less ob- 
scured case do not bear a strong resemblance to the AGN 
fraction curve. 

In the highly obscured case, the fis/fe slope is reddest 
when the AGN is most luminous, a somewhat counter- 
intuitive result since the 6/im continuum is generally be- 
lieved to be enhanced by AGN torus emission. This re- 
sult likely owes to the fact that the time when the AGN 
fraction is greatest corresponds to the time of greatest 
obscuration, when the 9.7/xm silicate absorption feature 
is deepest. This obscuration would lead to significant ab- 
sorption of the observed 6/im continuum, reddening the 
fib/ fe color. However, this may result from the simula- 
tions not having sufficiently clumpy ISM to allow for lines 
of sight directly to the AGN torus. This slope is also red- 
dened at t%, the time of peak SFR activity, as expected 
due to enhanced HH-region type emission. For the less 
obscured case, the /15//6 color is reddened throughout 
the merger {t\-t^, likely the result of the enhanced SFR, 
but lacks a sharper reddening at the time of coalescence. 
This supports our view that the "coalescence" reddening 
seen in the "highly obscured" case is a direct result of 
galaxy dust attenuation. 

Similarly, the /30//15 slope is effectively constant for 
the duration of the merger for the less obscured case, but 



is significantly reddened at the time of coalescence in the 
highly obscured case. At a given time, this slope is bluer 
when the AGN source is more luminous (e.g., AGNxlO), 
and in contrast to fxs/fe, /30//15 does not experience a 
sharp peak at £3 for the most powerful AGN in the highly 
obscured merger. 

Overall, our simulations indicate that the dependence 
of any of these mid-IR spectral diagnostics on the AGN 
fraction is time-dependent and non-linear and is affected 
by the properties of the host galaxy. The dependence on 
the overall level of obscuration, as parametrized by the 
silicate absorption depth, is to some degree stronger. For 
example, these indicators vary much less when there is 
less obscuration. 

When the highly obscured and less obscured simula- 
tions have the same AGN fraction, the highly obscured 
case shows much sharper coalescence-stage SED redden- 
ing and deepening of the silicate feature. These effects on 
the SED are therefore the direct result of the host galaxy 
dust attenuation. This supports ear l ier tentative obser- 
vational evidence (|Lacv et all l2007t iSaiina et ail 12007ft 
that in at least some AGN-dominated, deep silicate ab- 
sorption sources, the absorption could arise in the host 
galaxy rather than the AGN torus. 

4.2. The Uncertainty of Dust 

In this section, we address the effect of varying both 
the dust structure in the ISM (clumpy vs. smooth) and 
its composition (Milky Way- like vs. SMC-like). In Sec- 
tion 2.4, we described the two ISM treatments we apply 
to the highly obscured merger. Broadly speaking, the 
"default ISM" case assumes the dust and gas associated 
with both hot and cold ISM phases are smoothly dis- 
tributed. The "alternate ISM" assumes the cold phase 
gas is in clumps sufficiently dense to have a negligible 
volume filling factor, so that the dust it contains does 
not contribute to the overall attenuation. In Figure 21 
we explore these effects on the mid-IR SEDs. 

When switching to the Alternate ISM assumption, the 
magnitude of rg.7 is reduced by ~ 0.3 dex, and the color 
/30//15 by ~ 0.7 dex, while r 7 . 7 and fi 5 /f 6 are mostly 
unchanged. In addition, the total IR luminosity (first row 
of Figure U) is negligibly affected by this change. There- 
fore the differences from this ISM structure assumption 
arise only because the source energy is redistributed in a 
different way, and not to changes in the total amount of 
energy that is attenuated. 

With SMC or MW dust, the IR SEDs at a given 
time during the starburst have similar slopes /30 / /15 and 
/15//6 (see also, Figure [T]). Furthermore, Ljr, a mea- 
sure of the amount of attenuated light, is nearly iden- 
tical in these cases. By contrast, the SEDs show very 
different spectral characteristics r7.7 and rg.7, reflecting 
the difference between MW and SMC grain composition. 
The SMC dust is assumed to contain relatively more sil- 
icate than carbonaceous grains, leading to virtually non- 
existent PAH emission r7.7 and correspondingly stronger 
silicate absorption rg.7. 

4.3. Viewing Perspective 

In addition to the (evolving) source SEDs, the observed 
SED depends on the evolving relative distributions of 
sources and dust. This evolution changes the relative 
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Figure 4. We compare the dust models in the same simulations 
and quantities as Figure [3] For both bursts, we show our fiducial 
models in solid black curves, the SMC dust grain model in the 
dashed purple curves, and the fiducial models without galaxy dust 
in the light blue curves. Vertical lines highlight the four character- 
istic times t±-t4 from Figure [TJ corresponding to pre-burst, peak 
starburst, peak AGN, and post-burst. In the deeply obscured case, 
we also plot two alternate models: the SMC dust grain model and 
the clumpy ISM structure model. The powerful AGN signal in 
rg.7 and T7.7 at tz (Figure [3j can be replicated, at normal AGN 
strength, with the silicate-heavy SMC dust grain distribution. In 
several rows, the solid gray (patterned blue) boxes subtend the 25- 
75 % range in the distrib ution of local ULIRGs (quasars) compiled 
bv lVeilleux etaT] J2009D . 

obscuration of AGN and stellar sources, affecting the be- 
haviors of mid-IR indicators depending on whether the 
luminosity is dominated by star formation or AGN, be- 
cause the AGN is generally much more highly obscured 
than the stellar component. 

In Figure [5] we examine the time-dependence of vari- 
ation in the mid-IR quantities with viewing angle, fo- 
cusing on the highly obscured merger. The camera an- 
gle scatter is much lower (median S < 0.05 dex) for 
the less obscured merger. We show the same assump- 
tions regarding AGN strength and ISM models from Fig- 
ure HI and we plot the scatter S, in dex, of the mid-IR 
quantities as a function of time. We define S(q) to be 
the normalized median absolute deviation (MAD, e.g., 
IMosteller & Tukevl[l977l : IBeers et alJ[l990l ) of the quan- 
tity q among the seven Sunrise viewing directions at a 
particular time during the simulation. First we calculate 
the median M q of q, and then the seven values \q — M q \, 
the median of which we define to be the MAD of q. We 
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Figure 5. Scatter S, in dex, between viewing angles for the in- 
frared quantities of the highly obs cure d case from Figures \3\ and [4] 
S is defined in the text of Section I4.3I Vertical lines highlight the 
four characteristic times t\-t^ from Figure[JJ corresponding to pre- 
burst, peak starburst, peak AGN, and post-burst. The scatter S is 
largest when the AGN is strongest, ~ tz—tn, because the obscura- 
tion varies significantly among the different sightlines to the AGN. 
The camera angle scatter for the less obscured burst (not shown) 
is much smaller, with S always below 0.1 dex, and its median value 
below 0.05 dex. 



define S(q) = MAD/0.67 so that if q is drawn from a 
Normal distribution with standard deviation cr, then S is 
an unbiased estimator of a . 

At t < t 2 , S < 0.1 dex for all models, so the SED 
shape and spectral features are similar in all directions. 
Generally speaking, the scatter S reaches a maximum 
(~ 0.2-0.4 dex) for conventional AGN diagnostics dur- 
ing £3 to £4. These times are also when the AGN lu- 
minosity fraction and IR luminosity are highest. Thus 
the increased viewing-angle scatter implies a large AGN 
contribution to the IR SED, and that this contribution 
depends sensitively on the distribution of dust in the 
host galaxy along the line of sight. Interestingly, r7.7 
experiences much less viewing angle scatter than rg 7l 
S{r 7 . 7 ) ~ 0.0 versus S{r g . 7 ) ~ 0.2 dex. Turning off AGN 
emission leads to essentially no viewing angle variation 
in both quantities. 

Figure[S]presents the SED from seven viewing angles at 
time £3 (peak AGN) for our highly obscured merger, with 
AGN x I and fiducial ISM models. This gives us a clearer 
physical picture of what is driving changes to the mid- 
IR diagnostics during coalescence. The generally higher 
level of obscuration at these times leads to the enhanced 
silicate absorption, causing the dip in median rg. 7 at £3 
(Figure [3]). However, at a fixed time during this period 
of very high attenuation, like the one shown in Figure |H1 
the silicate absorption feature is stronger along lines of 
sight in which the AGN is less obscured. In other words, 
the silicate feature only correlates with AGN attenua- 
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Figure 6. Infrared SEDs of the AGN X 1 highly obscured burst 
simulation at time tz (AGN peak) observed from each of the seven 
viewing angles. The IR spectrum varies owing to the anisotropy of 
extremely optically thick material surrounding the central engine 
(starburst or AGN). Along a given direction, this dust absorbs IR 
light and reemits it into other lines of sight. 

tion until the level of attenuation is so extreme that the 
emergent stellar SED dominates the mid-IR emission as- 
sociated with the AGN. This galaxy dust also attenuates 
both the T7.7 emission and IR continuum. 

This ultra-obscured example is quite rare in the simula- 
tions, as we find in Section[5J However, this example also 
demonstrates why, at a given time in our highly obscured 
merger, the PAH features experience less viewing-angle 
dependence than the silicate absorption depth. Regard- 
less of the attenuation along a given line of sight, the 
PAH emission and mid-IR continuum are obscured by 
the same amount because they are emitted by the same 
dust. By contrast, the silicate absorption feature de- 
pends both on the amount of obscuring dust and also 
the level of the mid-IR continuum, which may itself be 
heavily obscured. 

4.4. Intrinsic AGN Mid-IR Emission 

In Figure [7] we consider the effect of intrinsic mid-IR 
AGN SED shape on the final observed SED. We demon- 
strate that the intrinsic AGN mid-IR spectrum has no 
impact on the observed mid-IR SED at £2-^3 of the 
highly obscured merger, during the starburst and AGN 
peak at final coalescence when the central densities are 
hig hest. We use our def ault AGN SED, the template 
by iHopkins et ail (l2007cD . as well as two templates by 
iNenkova et al.l (|2008|) . These templates reflect obscura- 
tion by clumpy torii inclined by i = 0° and i = 90°, 
spanning the range in mid-IR shape they derived. 

Dust obscures this central source at wavelengths as 
long as A ~ 20-50/xm. Thus the flux at A ~ 2 — 5/xm, 
conventionally attributed to AGN torii, only indicates di- 
rect AGN emission when the central source is sufficiently 
unobscured (e.g., t\ and £4). In the less obscured simu- 
lation, the AGN torus SED has a large effect throughout 
the merger (Figure [TJ. 

5. OBSERVATIONAL COMPARISON 

5.1. 2- Dimensional Diagnostics 

In Section 01 we presented the time evolution of the 
observed mid-IR indicators and noted that their depen- 
dence on Lagn/^boi is both time-dependent and non- 



Figure 7. We explore the sensitivity of the final infrared SED 
to the assumed intrinsic AGN (or AGN+torus) spectrum, focusing 
on the highly obscured simulation. Solid lines show the emergent 
SED from an example viewing direction with three assumptions 
about the intrinsic AGN SED. In this figure we use the manually 
boosted AGN (xlO) input SED, and so the effect sizes shown here 
should be considered upper limits. Dashed lines show the intrinsic 
stellar+AGN source emission before the dust RT calculation. Black 
lines use the mean templates from Richards et al. (2006a), while 
blue and yellow lines use i = 0° and i = 90° clumpy torus models 
from Ncnkova et al. (2008), respectively. The choice of intrinsic 
AGN SED matters least at tz, when the SMBH's contribution to 
the total luminosity is largest. The AGN fraction at times ti~ t4 
are all > 70% and yet the appearance of the mid-IR SED changes 
dramatically over these several 10 s yr. The influence of the intrinsic 
AGN spectrum on the final SED at 5fim increases from ~ dex 
at tz to ~ 1 dex at t±. 

linear. We showed that our mid-IR spectral diagnostics 
strongly depend on the dust model and can be highly 
scattered as a function of veiwing angle. How could we 
hope to observationally discern the level of AGN in a 
given galaxy based on its mid-IR spectrum, given that 
typically no information on the merger stage is available, 
the dust model is even less constrained, and we will never 
have more than one viewing perspective? 

In Figure [U we present mid-IR diagnostic diagrams 
that have been used in the literature for the purpose 
of disentangling the role of AGN and starburst ac- 
tivity in dusty galaxies. Specifically, using our def- 
initions of the PAH strength, silicate feature depth, 
and mid-IR colors, we present s lightly modified ver - 
sions of the "Lauren t" diagram (jLaurent et al.l [2000) , 
the "Spo on" diagram dSpoon et al.ll2004D . the "Veilleux" 
diagram dVeilleux et al.l 120091 ) . and the "Lutz" diagram 
(jLutz et al.l2004l ) . The last diagram presents the relation 
between the 6/im luminosity and X-ray luminosity for 
unobscured or mildly obscured AGN. The left-hand col- 
umn of Figure [8] presents the highly obscured simulation, 
the middle column presents the less obscured simulation, 
and the right-hand panel presents real galaxy data that 
we discuss in more detail in the following section. 

From the simulations, we plot the desired IR quantities 
into a grid of 25x25 squares, and calculate the median 
value of Lagn/^boi for each bin (including all possible 
viewing angles), which are color-coded as indicated by 
the colorbar. From this coding, we can identify which 
regions, if any, isolate powerful AGN in the models. 

5.2. Comparison Between Models and Data 

I Veilleux et al.l (|2009l ) anal yzed the standard ev olution- 
ary scenario for ULIRGs (| Sanders et al.l 119881) , which 
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Figure 8. We compare the two simulations to observed samples of ULIRGs at various redshifts. The four rows each present an example 
diagnostic diagram that has been proposed and examined in the literature. Within each row, the left plot shows the "highly obscured" 
model points colored according to the median AGN fraction £ a gn/£bol within each grid cell: see the text for a description of this procedure, 
and refer to the legend in the upper center of the figure for the translation. The middle panel of each row shows the same mapping for 
the "less obscured" simulation. The right panel in each row plots observed local V09 and high-z S12 ULIRGs. When the rg. 7 feature is 
saturated, we plot a b rown asterisk. In t he "Lutz diagram" , the gray shaded region shows the la range of intrinsic X-ray-to-mid-IR ratio 
for Type 1 AGN from Bauer ct al. (201(3). For the models here we assume MW dust. 



broadly speaking proceeds through three stages. In 
"Stage 1" , the IR emission of a gas-rich merger is domi- 
nated by star formation, leading to large PAH EWs and 
relatively small silicate absorption EW. As the merger 
proceeds to "Stage 2" , r 9 .7 increases as the central source 
is buried in the gas- rich starburst, and the PAH emission 
decreases relative to the IR continuum that is boosted 
by the AGN, which remains sub-dominant to the star- 
burst. Finally, at "Stage 3", the silicate absorption 
and PAH emission both subside as the starburst fades, 
the obscuring medium is cleared or consumed, and the 
UL IRG is comple t ely or nearly AGN dominated. Over- 
all. lVeilleux et all (|2009fl concluded that AGN contribute 
~ 35 - 40% to the IR luminosity of ULIRGs, and that 
ULIRGs experience significant scatter from this general 
pattern. 

In Figure [U the highly obscured model points roughly 



fill the space spanned by observed sources. Our pro- 
cedure does not mandate such consistency, and so we 
interpret this as one important check of the modeling 
technique. In particular, ignoring other diagnostics (such 
as SED slopes), this simulation reproduces the observed 
spread in rg.7 and r^.f. In Table [IJ we summarize a 
few numerical results of our deeply obscured simula- 
tion. We divide the merger into three stages, t\-ti 
(pre-coalescence), t 2 -<3 (coalescence), and t^-ti (post- 
coalescence), and report the median values of iagn/^bob 
rg.7, and r7.7 calculated in those periods. 

Generally, the SED properties along the model time- 
line follow those of the merger stages by V09, albeit with 
an expectedly high amount of scatter. Local ULIRGs 
have a large amount of silicate absoprtion (rg.7 ~ 0.7), 
suggesting they may correspond to the modeled points 
at ti < t < £3 where rg.7 peaks. However, the AGN frac- 
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Table 1 

Selected Diagnostics. Deeply Obscured Case 



time range a logr7.7 
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-0.37 


0.00 


ALiiN xi, iviw dust 






tl <t<t2 0.57 


-0.33 


0.06 


IZ < i < to U.4o 


—0.62 


0.33 


£3 < t < £4 0.39 


-0.43 


0.26 


a^jIn xi, oivio qusl 






tl < t < t2 0.12 


-0.44 


0.06 




— 1.09 


0.33 


t3<t<t4: -0.03 


-0.70 


0.26 


AGN xlO, MW dust 






tl < t < t2 0.40 


-0.37 


0.43 


t2 < t < t3 0.04 


-0.95 


0.86 


t3 < t < t4 0.30 


-0.40 


0.81 


Combined, MW dust 






tl <t <t2 0.57 


-0.32 


0.36 


t2 < t < t3 0.45 


-0.65 


0.70 


t3 < t < t4 0.39 


-0.38 


0.58 



a These time ranges are 0.154, 0.056, and 0.252 XlO 9 
yr long, respectively 
b Median values 



tions in models depend sensitively on the dust model and 
AGN strength choices, and so it is unclear how well they 
may correspond to the ~ 30% reported in V09. This 
value is consistent with our AGN x 1, MW dust simula- 
tion, but is much lower than the ~ 70% from the "com- 
bined" model set, in which each of the AGN strengths 
are given equal weight. The more realistic choice is un- 
certain, but the AGN x 1 case has trouble reproducing 
observed objects with rg.7 < — 1. However, SMC dust at 
lower AGN strengths can provide that feature. There- 
fore, it is very difficult to determine the path of the AGN 
fraction in these objects. 

In Figure HI swaths of strong and weak AGN exist, but 
most of these projections do not isolate AGN in simple 
ways. For example, while it is true that sources with 
strong starbursts and weak AGN typically have stronger 
PAH emission and redder /30//15 as discussed in Veilleux 
et al. (2009), strong AGN occupy all values of /30//15, 
/15//6J an d rg.7. Furthermore, the PAH strength de- 
pends sensitively on the dust model we assumed (see Fig- 
ureS]), and the r^n values of low £ agn /Lboi points depend 
on the obscuration level of the simulation. L agn /Lboi ~ 
squares exist at a wide range of SED slope, rg.7, and r-j.T, 
confounding potential clean estimates of L a gn/-^boi in 
our massive hyper-LIRG example (the "highly obscured" 
merger). These model points occur after £4, and after 
most of the SF and AGN activity has died away. How- 
ever, these sources are still bright enough in the IR to 
be considered ULIRGs, and so we have included them in 
this analysis. 

In addition, we compare the relative fluxes emerging at 
hard X-ray and mid-IR wavelengths in the bottom row 
of Figure jH We compute m odel X-ray fluxes as in Sec- 
tion H31 iBauer et all (|2010t ) found that z > 1 ULIRGs 
that are otherwise thought to be AGN-dominated (with 
low ry 7) are undetected at X-ray wavelengths, suggesting 



that they are at least mildly Compton-thick. Moreover, 
they presented Compton-thick examples from the litera- 
ture at z ~ 0, which li e ~ 2 dex below t he relation for 
less obscured AGN (c.f.. lHonig et al.H2010D . Many model 
galaxies with powerful AGN lie in the same region of this 
diagram, and are reflected in the dark blue squares mixed 
in with the lighter squares at Lq ~ 10 11 -10 12 Lq. These 
are the same sources in the low-r7.7 (high rg.7) tail of 
dark points in the upper panels, implying that this class 
of Compton-thick sources is a short-lived phase of our 
highly obscured model. 

There are some properties for which the model val- 
ues are not realized by observed sources. Generally 
speaking, the differences are most apparent in the mid- 
IR SED slopes /30//15 and f 1G /f 6 . The most AGN- 
dominant model points are pa 0.5 dex redder in /15//6 
than any observed object. In bulk, the model /30//15 
values are larger by w 0.3 dex than the observed distribu- 
tion. And there are few model galaxies with log rg.7 ~ 
and logr7.7 ~ 0, in contrast to those observed. Ho wever , 
this is sensitive to our dust grain model (Section 16.21) . 
and may also reflect our not having modeled AGN with 
no galaxy obscuration. Our calculations are necessarily 
uncertain, so in Section 16.31 below, we summarize some 
of these theoretical limitations and potential implications 
of these mismatches. 

As we have seen in Section |4l when sufficiently ob- 
scured, a powerful AGN has an S ED shape the same as 
a dusty pure starburst (see also I Younger et al.l [2009b; 
INaravanan et al.ll201 0b). Figure [5] shows high L agn /Lb \ 
squares occupying a large fraction of the ?"7.7, rg.7, 
/30//15, and /15//6 values. We also saw that the model 
SEDs can change drastically on timescales ~ 10 8 yr. A 
further challenge is that the implied £ agn /Lboi, given a 
set of SED properties, depends sensitively on assump- 
tions about the dust model and AGN strength. An ideal 
indicator is one that identifies powerful AGN regardless 
of the host galaxy's content or state. An example is hard 
X-ray emission which can penetrate even very large op- 
tical depths, and in Section [SJ we use our simulations 
to suggest other idealized indicators that have not been 
easily accessible with current data, but may be useful 
for future mid-IR studies with the James Webb Space 
Telescope {JWST). 

6. IMPLICATIONS 
6.1. An Ideal Indicator of L agn / Lbol 

Ultimately, we seek an indicator that can be used to 
determine the AGN contribution to IR-luminous galaxies 
across the full range of obscuration by galaxy dust. Near- 
infrared colors effectively identify AGN when their IR 
SED resembles a power-law, leading to enhanced emis- 
sion beyond the stellar bump, e.g., > 2/im, and a high 
value of our example near-IR color f4.5ff1.e- This tech- 
nique has been used to successfully identify high-redshift 
AGN-dominated ULIRGs using Wide- field Infrared Sur- 
vey Explorer (WISE) photometry (e g IStern et al1l2012t 
lEisenhardt et alj|2012t IWu et aLlfeoia . 

This near-IR reddening operates in powerful AGN for 
a wide variety of attenuation. However, as we saw in 
Figure [71 an extremely high dust column density can, on 
short timescales, change the near-IR SED shape from 
a power-law AGN to one that resembles a starburst- 
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Figure 9. Following Figure[8] we plot panels motivating the diag- 
nostic D: see text for definitions of Dmw an d -DgMC- The left two 
columns are our simulated points, colored by L agn /Lt,ol! f° r our 
two representative levels of obscuration. The top row assumes MW 
dust, and the bottom assumes SMC dust. In contrast to the other 
quantities from Figure [8] diagnostic D straightforwardly estimates 
^agn/^bol f° r both low and high optical depth, assuming we know 
the dust model. We note that diagnostic D requires information 
at rest-frame wavelengths 1— lOfim, rendering it somewhat imprac- 
tical to measure using any single current instrument. However, 
diagnostic D will be accessible to the James Webb Space Telescope 
for sources at z < 1. 

dominated ULIRG. In this situation, the information 
about the dust attenuation contained in rg.7, a fairly ob- 
vious feature when present, may be helpful. Consider a 
combination of rg.7 with near-IR color /4.5//1.6, such as 



Dmw = \ (log r 9 . 7 + 0.32) 



logy— 

jri.e 



0.5 



The effectiveness of combination Dmw is evident in the 
top row of Figure |H1 the map between color/shading 
(iagn/iboi) and these axes, /4.5//1.6 and rg.7, appears 
to be a well defined and slowly varying function of few 
parameters. 



'.2. Constraining the Dust Properties 
We note that this depends on our dust model: 



the 



bottom row of Figure [9] shows the same diagnostic as- 
suming SMC dust instead of MW dust. Importantly, a 
nice mapping in the above sense still exists, it is just 
skewed toward higher values of rg.7. We define an alter- 
nate diagnostic, 



D 



SMC 



^/(logrg.7 + 0.32) 2 



logy- 
Jl.6 



0.5 



centered around the same point as Dmwj but with a 
different axis ratio and scaling. 

In order to select between these estimators, we must 
first constrain the dust properties; in other words, we 
must jointly constrain the dust and i a gn/^boi- One op- 
tion is to use the PAH emission: assuming SMC dust, 
r 7 7 ~ for all model points instead of ~ 0.3 for MW 
dust (see Figure 2]). We demonstrate this procedure in 
Figure [TUJ while Dmw does not yield a low-scatter pre- 
dictor of £ agn /Lboi when applied to a system with SMC 
dust, it does broadly select stronger AGN. Moreover, at 



median log L a <, n /L bol 



□ 



-1.3 -1.0 -0.7 -0.4 -0.1 




MW 



Figure 10. PAH strength versus our combination diagnostic 
Amy, which contains information about obscuration rg.7 and 
near-IR color, providing a rough estimate of L agn /Lb i at both 
low and high mid-IR optical depths. If we have incorrectly assumed 
MW dust when the appropriate model is more like the SMC dust, 
then Figure l9l implies that the Dmw - ^agn/Aiol relation will be 
skewed and have a large scatter. However, at fixed high DmWi the 
PAH strength can roughly select the appropriate dust model, and 
therefore the appropriate predictor of L a gn/ £*bol : fMWi DsmCi 
or an interpolation. 

D > 0.2, sources with these different dust models sepa- 
rate based on their PAH strength, 7-7.7, permitting either 
a direct choice between the formulae for Dmw and Dsmc 
or an interpolation between them. This suggests the fol- 
lowing procedure: 

1. Estimate Dmw- 



2. Estimate the PAH strength (here, tt.t) at D 



MW- 



3. Select dust model X based on the value obtained 
in step 2. 

4. Calculate Dx and use it to estimate £ a gn/-^bol- 

With a known dust model, diagnostic D robustly pre- 
dicts Lagn/iboi regardless of the physical conditions of 
the ULIRC source. We compare the ability of diagnos- 
tics Dmw and Dsmc to predict £ agn /Lt>oi to that of hard 
X-ray fluxes in Figure QT] 

Using a least absolute deviation method, we fit the re- 
lation iagn/^bot = a+bQ for three quantities Q and each 
of our two ISM models for the highly obscured simula- 
tion, and compile the results in Table [2l To measure the 
scatter in the correlations' residuals, we used the same 
scale estimator S defined in Section [4~3l which is gener- 
ally insensitive to outliers. In both cases, La— 10 keV and 
each D predict L agn /Lboi with roughly the same median 
scatter, S ~ 10%. However, we find that L2-10 kcV ex- 
periences a failure rate (> 3a) two to three times higher 
than D in the highly obscured calculation. Both quan- 
tities predict L agn /£boi equally well under the alternate 
ISM treatment and for the less obscured candidate (not 
shown), with S < 10% and failure rates < 1%. 

Several simulated points at L agn /iboi ~ 0.8 fail in both 
L2-10 kcV and Dmw- In these cases, corresponding to 
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Figure 11. We show how diagnostic D correlates with iagn/^bol m our simulations in the right three panels. We show here only the 
highly obscured merger with our default ISM model, but the trends are similar for the alternate ISM model and the less obscured merger, 
as summarized in Table [2] As a comparison, we plot the hard X-ray flux through the 2 — 10 keV band in the left panel. If the correct dust 
model is used, D accurately predicts iagn/^bol as well as Z.2— 10 koV- In th e rightmost panel, we apply -Dmw when the sou rce has SMC 
dust instead, demonstrating that the scatter increases by ~ 70%. This issue can be mitigated by using, for example, Figure [T0l to choose 
the appropriate indicator D. 

A concern with the present method is one of "progeni- 
tor uncertainty" , whereby we cannot be sure how our few 
idealized simulations fit into a cosmological context. We 
have mitigated this issue by presenting two rather dif- 
ferent ULIRG candidates, by varying the relative power 
between the AGN and star formation, and by varying 
the assumptions acting on dust obscuration and emis- 
sion. Over time, we seek to build a library of model 
SEDs for a much wider variety of systems, enabling more 
robust studies across IR-luminous galaxy populations in 
the future. Moreover, we have focused here on AGN sig- 
natures at mid-infrared wavelengths, and infor mation at 
longer wavelengths (c.f.. lYounger et aT1l2009aft . such as 
the location of the far-IR peak, could further constrain 
the AGN fractions in such sources (e.g., S12). 

The simplest interpretation of Figure [8] is that ana- 
logues of our "highly obscured" and "less obscured" ex- 
amples are both realized in observed sources, in addi- 
tion to both MW and SMC dust: all four models are re- 
quired to cover the wide observed distributions of mid-IR 
color, PAH emission, and silicate absorption. A higher 
ratio of silicate to carbonaceous dust grains increases 
rg.7 while shrinking rj^ (Figures |4|) , without changing 
/4.5//1.6 (Figure O, permitting estimates of L &sn /L ho \ 
under a variety of physical conditions, dust models, and 
optical depths. Although different ISM or dust models 
can produce broadly similar mid-IR trends, a separate 
measure of the AGN power, such as X-ray flux, is essen- 
tial to further constrain these models. Furthermore, the 
recent discover y of high-redshift , hyp e r-luminous, hot- 
dust ULIRGs (|Eisenhardt e"t~ai1 l20ia IWu et all I20H 
suggests that we may lack a key phase in these few sim- 
ulations, a transition between SEDs represented by our 
"highly obscured" and "less obscured" examples. 

Several assumptions regarding the dust radiative trans- 
fer may lead to discrepancies between our model sources 
and real sources. In particular, from Figure [8] we have 
identified two failures to which these assumptions likely 
contribute: 1) the mid-IR colors of the most AGN- 
dominated highly obscured sources are redder than their 
real counterparts (although bluer at A < 5/xm; see Fig- 
ure^, and 2) we do not simulate enough sources with lit- 
tle signs of PAH emission or silicate absorption: rg.7 ~ 
and r7.7 ~ 0. 

The IR calculation in Sunrise distributes the energy 
absorbed by small PAH grains equally between emis- 



Figure 12. Synthetic photometric observations by MIRI on 
JWST, demonstrating the ability of future mid-IR technology to 
measure AGN power in systems with a wide range of obscuration. 
Here we assumed MW dust. 

the most highly obscured lines of sight in our AGN x 
10 calculation from Figure [51 the mid-IR is so heavily 
attenuated that no rg.7 signature from the central source 
survives in the observed SED (see also Section |4~31 . 

Diagnostic D requires information about the galaxy's 
stellar mass (rest-frame A ~ 1-2/xto), limiting its imme- 
diate application to large new samples. However, the 
rest-frame wavelength range D requires will be covered 
by the James Webb Space Telescope {JWST), an ob- 
servatory whose M id-Infrared Instrument (MIRI) (e.g., 
I Wright et al. 2004) permits measurements of T7.7 and rg.7 
for sources at z < 1. The tight correlation we find be- 
tween D and I/ a gn/iboi reflects the possible benefit in 
combining information about rg.7 and r7.7 with near-IR 
color diagnotics from JWST itself (e.g., iMessias et al.l 
|2012| ). or by matching to sources from other IR sur- 
veys (e.g., iStern et all 120051 I2012T ). For the two ex- 
ample SEDs in Figure [121 whose known AGN fractions 
are both L a gn/^boi = 0.33, this procedure estimates 
iagn/^boi = 0.33 ± 0.09 and L asa /L hol = 0.38 ± 0.09 
(uncertainties are minimum theoretical scatter) for the 
high and low obscuration examples, respectively. 

6.3. Limitations and Improvements 
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Table 2 

Ability to Predict iagn/^bol 



Quantity, Q a 


scatter, S h 


intercept, a 


slope, b 


Failure % c 


Default ISM 










logL 2 -10 kev/^O 


0.12 


-4.51 


0.47 


6.1 


D MW , MW dust 


0.09 


-0.04 


1.07 


2.6 


D S MC, SMC dust 


0.12 


-0.01 


0.93 


0.5 


D MW , SMC dust 


0.14 


0.05 


0.72 


0.8 


Alternate ISM 










logL 2 -io kev/^0 


0.09 


-4.59 


0.47 


0.8 


Dmw, MW dust 


0.10 


-0.04 


0.96 


0.9 


-Dsmc, SMC dust 


0.09 


0.00 


0.86 


0.0 


Dmw. SMC dust 


0.10 


0.00 


0.98 


1.3 



a L am /L hol = a + b(Q) 

b S = MAD/0.67, see Section [431 

c Percentage of points with L agn /L bo i further than 3cr away from the 
fitted relation. In practice, for these cases the indicator predicts an 
iagn/^bol smaller than the true value owing to extreme attenuation. 

sion in thermal equilibrium and in modes giving rise t o 
the distinctive mid-IR PAH lines (jJonsson et al.ll2~01Q[) . 
Therefore, any change to this assumption or more sophis- 
ticated calcu lation, such as thermally fluct uating small 
grains (e.g., iGuhathakurta fc Drainel fl989l) . will affect 
the dust continuum shape and PAH strength. In ad- 
dition, we do not currently model dust destruction, and 
our model for dust production is crude, both key areas 
where near-future progress may be tenable. By not al- 
lowing our radiation to destroy dust grains or alter the 
dust grain distribution (or changing more generally the 
gas distribution owing to realistic feedback), we may be 
preventing ourselves from obtaining truly "power-law" 
AGN SEDs with no mid-IR dust features that may arise 
after AGN feedback disrupts the ISM. 

We have purposefully focused on how galaxy dust can 
affect the fixed SMBH source SED (with minor variations 
in Section [O]) , rather than how the physics on smaller 
scales play out in observations of ULIRGs. While high 
levels of galaxy dust can reproduce a number of prop- 
erties of obscured quasars, the evolution of the central 
regions during these times of maximum obscuration will 
depend on the behavior of gas structures very close to 
the SMBH. This behavior may set the minimum scat- 
ter to be observed in correlations between these i ndica - 
tors. Using multiscale simulations, Hop kins et al.l (|2012[ ) 
showed how gravitational instabilities propagate down 
to scales of ~ 1-100 pc during large-scale gas inflows. 
They demonstrated that a number of properties of ob- 
scured AGN attributed to a torus are explained by such 
dynamically-evolving accretion structures. Therefore we 
are motivated to develop methods of integrating knowl- 
edge of the impact of these different scales on the emer- 
gent SEDs as a means for further interpreting the evolu- 
tion of IR-luminous galaxies. 

7. SUMMARY AND CONCLUSIONS 

We have simulated two representative starbursts meant 
to bracket the plausible conditions in observed ULIRGs. 
We studied the variation of the observed mid-IR SEDs 
as a function of viewing angle, source SED, dust grain 
distribution, galaxy ISM dumpiness, and AGN power. 
By calculating commonly used SED diagnostics along- 



side the modeled SMBH accretion rate, we studied the 
ability of these signatures to predict the AGN UV-IR 
luminosity fraction, Lagn/^boi, and developed a mid-IR 
estimator that can be measured for new samples with 
next-generation instruments such as JWST. Our key con- 
clusions are: 

1. Our modeling technique broadly reproduces the ob- 
served span of common mid-IR diagnostics used to 
disentangle AGN and starburst activity. 

2. Although generally consistent with previous in- 
terpretations, none of the indicators straightfor- 
wardly predict L agn /Z/boi because they depend non- 
linearly on the galaxy's properties, which can vary 
rapidly on timescales ~ 10 8 yr. 

3. Highly obscured sources can be more AGN- 
dominated than realized, even when the shape of 
their SED resembles a pure starburst. A large 
enough column density reprocesses all of the direct 
AGN light into longer wavelengths. 

4. The mid-IR SEDs of sources with a powerful, 
highly obscured AGN depend strongly on the di- 
rection from which the source is observed. Along 
some lines of sight, dust self- absorption shrinks the 
IR luminosity and obscures the mid-IR signatures 
of the buried AGN. 

5. Sources with a deep 9.7/zm silicate absorption fea- 
ture can be reproduced by models either with an 
extremely luminous AGN and MW galaxy dust, or 
with an average-strength AGN and SMC galaxy 
dust. 

6. It is possible to construct an idealized mid-IR in- 
dicator utilizing future JWST data to cleanly esti- 
mate iagn/iboi while simultaneously constraining 
the dust grain model for any level of obscuration. 

7. Remaining mismatches between the observed mid- 
IR color distribution and our models imply a need 
for calculating dynamically coupled dust grain pop- 
ulations and estimating the effect of obscuration on 
scales nearer to the SMBHs than we have consid- 
ered. 



This research has made use of NASA's Astrophysics 
Data System. The computations in this paper were run 
on the Odyssey cluster supported by the FAS Science 
Division Research Computing Group at Harvard Univer- 
sity. 
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